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ABSTRACT 

This paper presents a numerical study over a wide parameter space of the likelihood of the dynamical 
bar-mode instability in differentially rotating magnetized neutron stars. The innovative aspect of this 
study is the incorporation of magnetic fields in such a context, which have thus far been neglected in 
the purely hydrodynamical simulations available in the literature. The investigation uses the Cosmos++ 
code which allows us to perform three-dimensional simulations on a cylindrical grid at high resolution. 
A sample of Newtonian magneto-hydrodynamical simulations starting from a set of models previously 
analyzed by other authors without magnetic fields has been performed, providing estimates of the 
effects of magnetic fields on the dynamical bar-mode deformation of rotating neutron stars. Overall, 
our results suggest that the effect of magnetic fields is not likely to be very significant in realistic 
configurations. Only in the most extreme cases are the magnetic fields able to suppress growth of the 
bar mode. 

Subject headings: gravitation — hydrodynamics — instabilities — stars: neutron — stars: rotation 



1. INTRODUCTION 

Rotating neutron stars formed following the gravi- 
tational collapse of a massive stellar iron core or the 
accretion-induced collapse of a white dwarf can be sub- 
ject to various nonaxisymmetric instabilities depending 
on the amount and degree of differential rotation. The 
prospects of detection of gravitational radiation from 
newly born rapidly rotating neutron stars by the current 
and future generations of kHz-frequency, ground-based 
gravitational wave interferometers highly motivate the 
investigation of such instabilities. In particular, if the 
rotation rate is high enough and shows a high degree of 
differentiation, the star is subject to the so-called m = 2, 
dynamical bar-mode instability driven by hydrodynam- 
ics and gravity, with m being the order of the azimuthal 
nonaxisymmetric fluid mode e ±lmip . On the other hand, 
at lower rotation rates gravitational radiation and viscos- 
ity can drive a star secularly unstable against bar-mode 
deformation. These two flavors of the bar-mode instabil- 
ity set in when the ratio (3 — T/\W\ of rotational kinetic 
energy T to gravitational potential energy W exceeds 
a critical value (3 C . Early studies with incompressible 
MacLaurin spheroids in Newtonian gravity showed that 
the onset of the instability arises when (5 C ~ 0.27 and 0.14 



for the dynamic al and secular cases, respectively (Chan- 
drasekhar 1969). 



Improvements on these simplified analytic models have 
been achieved through numerical work. Newtonian and 
general relativistic analyses of the dynamical bar-mode 
instability are available in the literature, using both sim- 
plified models based on equilibrium stellar configura- 
tions perturbed with suitable eigenfunctions, and more 
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involved models for the core collapse scenario. Newto- 
nian hydrodynamical simulations have shown that the 
value of (3 C is quite independent of the stiffness of the 
equation of state, provi ded the star i s not s t rongly dif- 
ferentially rotatin g (see Houser et al. (19941; New et al. 
( |2000[ ); |Liu| ( |2002[ ) and referen ces therein ) . "On the other 
hand, relativistic simulations ( Shibata et al.|[2000[ ) have 
yielded a slightly smaller value of the dynamical insta- 
bility parameter (f3 c ~ 0.24 — 0.25). They have also 
shown that the dynamics of the process closely resem- 
bles that found in Newtonian theory, that is, unstable 
models with large enough f3 develop spiral arms follow- 
ing the formation of bars, ejecting mass and redistribut- 
ing th e angular momentu m. Further relativistic simula- 
tions ( |Baiotti et al.|[2007 1 have shown the appearance of 
nonlinear mode-coupling which can limit, and even sup- 
press, the persistence of the bar-mode deformation. It is 
also worth mentioning that, as the degree of differential 
rotation bec omes higher and more ext reme, Newtonian 



simulations (Shibata et al. 2002, 2003) have also shown 



that rotating stars are dynamically unstable against bar- 
mode deformation even for values of f3 of order 0.01. 

Whether the requirements for the development of the 
instability inferred from numerical simulations are met 
by the collapse progenitors remains unclear. Observa- 
tions of surface velocities imply that a large fraction of 
progenitor cores are rapidly rotating. However, it has 
been shown that magnetic torques can spin down the 
core of the pro genitor, leading to slowl y rotating neutron 
stars at birth ( |Spruit fc Phinney|1998[ ). The most recent 
computations ot the evolution of massive stars, which 
include angular momentum redistribution by magnetic 
torques and spin estimates of neutron stars at birth, lead 
to core collapse progenitors which do not seem to rotate 
fast enough to guarantee the un ambiguous growth of the 



et al 



canonical b ar-mode instability ( |Heger et al 
. [2006 1 . These estimates are m agreemen' 
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served periods of young neutron stars. However, rapidly- 
rotating cores might be produced by an appropriate mix- 
ture of high progenitor mass and low metallicity. Recent 
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estimates suggest that about 1% of all stars with masses 
larger than 10 solar masses wi ll produce rapidly rotating 
cores ( Woosley & Hegcr 2006 ) . Recent relativistic simu- 
lations loTlOargirsamj^^ mod- 



els carried out by |Dimmelmeier et al. ( 2008| ), which in- 
clude a state-of-the-art treatment ot the relevant physics 
of the collapse p hase and realistic precollapse rotational 
profiles (see also Ott et al. 20071, have shown that the 
critical threshold of the dynamical bar-mode instability 
is never surpassed, at least early after core bou nce. How- 
ever, a lar ge set of the models invest igated by |Ott et ah] 
(2007) and Dimmelmeier et al. ( 2008| ) do show that mod- 
els with sufficiently differential and rapid rotation are 
subject to the low—/? instability. In addition it is also 



worth mentioning the simulations of Baiotti et al. (2008) 



which show that hypermassive differentially rotating neu- 
tron stars form following the merger of low-mass binary 
neutron stars, when modeled as polytropes. The result- 
ing hypermassive star undergoes a persistent phase of 
bar-mode oscillations, emitting large amounts of gravi- 
tational radiation prior to its delayed collapse to a black 
hole. 

An important piece of physics that the existing nu- 
merical work has so far neglected is the presence of mag- 
netic fields. The amplification and emergence of strong 
magnetic fields in neutron stars from initially weak mag- 
netic field configurations in the pre-collap se stellar cores 



are currently under investigatio n (see e.g., Burrows et al. 
2007 Cerda-Duran et aL]|2008[ and references therein_ 



We note, however, that the weakest point of all existing 
magneto-rotational core collapse simulations to date is 
the fact that both the strength and distribution of the 
initial magnetic field in the core are unknown. 

In this paper, we show results from a detailed numer- 
ical study of the effects that magnetic fields may have 
on the dynamical bar-mode instability in differentially 
rotating magnetized neutron stars. In particular, we in- 
vestigate how sensitive the onset and development of the 
instability is to the presence of magnetic fields, as well 
as the role played by the magnetorotational instability 
(MRI) and magnetic braking mechanisms to alter the 
angular momentum distribution in the star and possi- 
bly suppress the bar-mode instability. Our study is mo- 
tivated by the potential astrophysical implications that 
the presence of strong magnetic fields may have for post- 
bounce core collapse dynamics and, in turn, for gravita- 
tional wave astronomy. 

The uncertainty of the strength and distribution of 
magnetic fields in collapse progenitors is reflected in 
our somewhat ad hoc parameterization of the field con- 
figuration through a large sample of equilibrium mod- 
els of rapidly and highly differentially rotating neutron 
stars. Our sample of Newtonian magnetohydrodynami- 
cal (MHD) simulations is based upon the set o f purely hy- 
drody namical models previously analyzed by |New et aL] 
( 2000[ ). The simulations are performed using the co yari 



ant (and adaptive mesh refinement) code Cosmos++ ( An- 
Fragilc 2003] |Anninos et"aT1p003| [20051 fragile 



mnos 



et al.|2005p which allows us to perform three-dimensional 



simulations on a logarithmically scaled cylindrical grid 
at high resolution. A number of equilibrium models of 
rapidly rotating stars with different values of the ro- 
tational instability parameter {(3 = T/lM^I) and mag- 
netic plasma beta {(3b = P/Pb), and different poly- 



tropic equations of state are constructed, introducing 
also different configurations for the magnetic field dis- 
tribution (of both poloidal and toroidal varieties) and 
field strengths. The equilibrium models are perturbed 
by seeding small random perturbations in order to initi- 
ate the onset of the bar- mode deformation. 

The organization of the paper is as follows. Section [2] 
discusses our basic formalism, numerical methods, diag- 
nostics, and the construction of initial data; Section K3] 
presents our results in two subsections, one for initially 
toroidal field configurations and one for poloidal. We 
conclude with a summary and discussion of our results 
in Section |U 

2. FORMALISM, INITIAL DATA, AND DIAGNOSTICS 

2.1. Cosmos++ Framework 

Although the Cosmos++ code now includes options to 
solve the full radiative and conductive MHD equations 
for Newtonian systems plus multi-species chemical and 
nuclear reaction networks, we only include here the sub- 
set of those equations necessary f or the current work. 
However, see Anninos et al. (20031 for a more complete 
de scription of the J Newto nian options including radiation, 
or jAnninos et al.j ( |2005[ ) for the general relativistic for- 
mulation and a discussion of the various energy formula- 
tion options available in the code. In the present study, 
we use the internal energy (and artificial viscosity, A.V.) 
formulation due to the robustness of the method in track- 
ing adiabats across thin kinematically dominated stellar 
atmospheres. The relevant equations are sufficiently dif- 
ferent than those published in our previous papers (sim- 
plified to Newtonian form and generalized to covariant 
curvilinear grids) that we write them out here for conve- 
nience: 



d(V9 P) 
dt 
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+ d t (JgS k V i ) = 
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di(Jg B l B k ) 



dt 

d{^gB k ) 



-V9d l [(P + P B )Sl + Ql] 
-yfgp d k <f> , (2) 

+ 9 s: (V5 e V) = -{PS} + (A) di{^/g V>) ,(3) 



Ot 



+ di{^g B k V 1 ) = ^gB l 8^ - g lk , (4) 



d t {Vd 9 ij d j $)=4nG^gp 



(5) 



where di = d/dt; 1 represents covariant derivatives in gen- 
eralized coordinates , and ^fg is the determinant of the 
spatial 3-metric g^ defining the coordinate system (cylin- 
drical for this work, with a logarithmically scaled radius 
to achieve greater resolution in the star's interior). Also, 
p is the fluid density, V k is the contravariant fluid ve- 
locity, Sk — pVk is the covariant momentum, e is the 
fluid internal energy density, Qj is the artificial viscos- 
ity tensor (here we use the u nsplit finite volume v ersion 
of the scalar viscosity from (Anninos et al. 20051 with 
k q = 2.0 and ki = 0.1), P is the fluid pressure, B k is 
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the contravariant magnetic field vector, Pb = B l Bi/8Ti 
is the magnetic pressure. <& is the gravitational poten- 
tial, which is found by solving equation (|5| with multi- 
pole boundary conditions that include up to 15 spherical- 
polar harmonics and all corresponding azimuthal mo- 
ments. In this work, we assume an ideal gas equation 
of state in the form P = (T - l)e. The MHD equa- 
tions are derived with the standard assumptions relevant 
for many astrophysical problems: the system is nonrcl- 
ativistic and fully ionized, the displacement currents in 
Maxwell's equations are neglected, the net electric charge 
is small, and the characteristic length scales are large 
compared to particle gyroradii scales. 

The scalar potential ip in the magnetic induction equa- 
tion Q is introduced as a divergence cleanser to maintain 
a divergence-free magnetic field (di(*Jg B l ) = 0). Op- 
tions are included in Cosmos to so lve any one of the fo l- 
lowing constraint equations for tp ( |Dedner et al.||2002 l: 



v 2 v = 



Of 



At 



-cldi{^gB% 



(G) 
(7) 
(8) 



which correspond, respectively, to elliptic, parabolic, and 
mixed hyperbolic and parabolic constraints. Here c p and 
Ch are user-specified constants used to regulate the fil- 
tering process and weight the relative significance of the 
hyperbolic and parabolic components. Also the time 
derivative in equation (j6| is approximated as a finite dif- 
ference with zero divergence at the initial time. For all 
of the calculations presented in this paper, we use the 
mixed hyperbolic and parabolic form with parameters 
Ch = 0.2c c flAx m i n /A£ and c 2 = 0.3c?,,, where c c a = 0.4 is 
the Courant coefficient, Aa; m i n is the minimum covariant 
zone length, and At is the evolution time step. 

Equations Q ~ ([5]) and (J8j) are solved in a modi- 
fied cylindrical coordinate system £ l = (j],z,(f>), where 
?7 = ln(ti7 +1) is a logarithmic radial coordinate used 
to concentrate resolution toward the interior of the star, 
with w = rsinO being the usual cylindrical radius. With 
this coordinate choice, the line element for the metric gij 
becomes 



dl z 



,2'/ 



drf + dz z + (e" - iy dtf , (9) 

and ^fg — e r, (e n — 1). We consider two different mesh 
resolutions, 64 3 and 96 3 , to address the robustness and 
relative convergence of our results. These are about opti- 
mal resolutions for three-dimensional simulations where 
a large number of parameters are to be explored, par- 
ticularly with active magnetic fields which increase sub- 
stantially the computational workload and suffer greater 
wave speed restrictions on the time step. 

2.2. Initial Data 

2.2.1. Rotating Polytropes 

We begin by constructing equilibrium models of 
rapidly rotating polytropi c stars using H achisu's self- 
consistent field technique (Hachisu 19861. For an ini- 
tially axisymmetric configuration with angular velocity 
ft = fl(zu) that depends only on the distance of the fluid 



from the rotation axis {w = rsin#), the equilibrium con- 
figuration satisfies 



(10) 



where ho and Co are constants, H = J p~ 1 dP is the 
fluid enthalpy, <f> is the gravitational potential obtained 
by solving the Poisson equation p|, and 



*{vj) = -i / \l 2 {zu)wdvu 



(11) 



describes the rotational profile of the star. For a poly- 



tropic gas with P = up = np 
H ={N +l)np 1/N 



(1 + 1/ N) wehave 



1/N 
v rmax,0 / 



1/N 



(12) 



where p max ,o is the initial maximum density. 

A variety of rotation profiles can be considered, such 
as rigid rotation (Q — const.), constant linear velocity 
(Q = Vq/xu), or constant specific angular momentum 
(0 = jo/vu 2 , where j is the specific angular momentum 
of the fluid). In our case, we choose a Maclaurin spheroid 
profile to com pare with the earl ier work of unmagnetized 
neutron stars (New et al. 1120001) 



Cl(zo) = ho 



1 



1 



m{w) 
M 



2/3* 



(13) 



where M equals the total mass of the star (spheroid), 
m(tn) is the mass interior to w, and 



h 



5J 
2M 



(14) 



can now be specified in terms of the total angular mo- 
mentum J of the star. 



The constants ho and Co in equation ( 10 1 are set by an 
appropriate choice of boundary conditions. Specifically, 
we require that p, P, and H vanish at the surface of 
the star. Any two points on the surface of the star can 
then be used to specify a solution; we choose a point on 
the equator (point ^4) and one of the poles (point B) by 
specifying the equatorial surface radius we and axis ratio 
zp/vje, where zp is the polar radius. The constants are 
then given as 



hl = - 



V(A) - *(S) 



(15) 



and 



Co = + h 2 Q * (A) 



(16) 

Finally we must specify T (or N) and n (or p maXi o) to 
close the system of equations. The solution proceeds by 
pressing an initial distribution for p; solving equations 
5|) and ( |ll| for $ and *, respectively; using equa tion 
TO I to set H and H max ; and then using equation (12 1 
to determine a new density distribution. This procedure 
is repeated iteratively until the solution converges suffi- 
ciently (for us, once Aho/ho, ACo/d, and AH/H are 
all < 10" 4 ). 

The initial data are solved in two-dimensional loga- 
rithmic cylindrical coordinates using the same polar and 
radial resolutions as the full three-dimensional grid that 
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we use for the production simulations, which we then 
map (i.e., rotate) azimuthally onto the three-dimensional 
grid. The quality of solutions is verified by computing the 
virial error (V.E.) as 



V.E. = 2T + W + 3 / PdV 



(17) 



where T is the total kinetic energy, W is the total gravi- 
tational potential energy, P is the thermal pressure and 
dV is the volume element. Virial errors in our calcula- 
tions are small: 4 x 10~ 4 for the T = 5/3 cases, 8 x 
for r = 2, and 1.6 x 1(T 3 for T = 3. Normalized 
to the absolute value of the potential energy, we find 
|V.E./W| < 3 x 10" 3 for all cases. 

Once the equilibrium density distribution peq is de- 
termined, we apply a small random perturbation of the 
form p{R, z, (ft) = Peq[1 + a>of(R, z, (ft)], where ao = 1CP 2 
is the perturbation amplitude and / is a random num- 
ber between -1 and 1. This perturbation serves as a seed 
for the bar-mode and magneto-rotational instabilities to 
grow. 

Because MHD codes have difficulty treating pure vacu- 
ums, it is necessary to initially fill the regions of the grid 
outside the star with a low density (/9fl 00 r = 10 -6 Pmax,o)> 
low energy (e floor = KpQ OOT / (r~l)), static (V 1 = 0) back- 
ground. During the evolution, any time the fluid density 
or energy attempt to drop below their respective floor 
values they are reset. 

2.2.2. Magnetic Fields 

In this work, we choose two idealized initial configura- 
tions for the magnetic fields, one purely toroidal and one 
purely poloidal. In all cases the fields are chosen to be ini- 
tially weak, so that in some sense they simply represent 
additional perturbations away from the initial equilib- 
rium state. However, it is now well understood that dif- 
ferentially rotating fluids with weak magnetic fie lds are 
susceptible to the MRI ( |Balbus fc Hawley|[l99l| . The 
only criterion required to trigger the MRI is 



dn 2 

dlntu 



< 



(18) 



a condition which is clearly met inside our model star. 
Thus, our initially weak magnetic fields can poten- 
tially become dynamically important during the evolu- 
tion (provided we have sufficient resolution to capture 
the MRI). 

The toroidal field case begins with B v = B z = and 



b 



0, Cj 



-A£ 2 /(2<r 2 ) 



-0.5 
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(19) 



where A£ = |£* — £i 00 pl ^ s the position offset relative to 
the location of the central strongest field loop £i oop = 
(??ioop) 0, <ft). The characteristic half-width of the loop is 
defined as a = we/8, and ?7i 00 p is set to either we/2 or 
we /4 so that the center of the toroidal loop is initially 
located at radius ri oop ~ 0.65zue or n oop w 0.28n7£. The 
normalization constant Ct is chosen to satisfy our partic- 
ular choice for /?B,min, the initially lowest value for P/Pb 
inside the star. The extra factor of e~ 0,5 is included to 



ensure that the magnetic field is initially confined to the 
interior of the star. 

The poloidal field case begins with = 0. The other 
two magnetic field components are specified from a mag- 
netic vector potential using the following divergence-free 
construction: 



>f3& = 







( 2 °) 

For the vector potential we use the solution correspond- 
ing to a current loop in the equatorial plane of radius 
n oop ( which we typic ally set to we 12) centered on the 
origin ( |Jackson||1975| ) 



A„ 



OO 

n=0 



(-l) n (2n- l)!!r< 
2 n (n + 1)! 7\ 



-2 P 2n+l(cOs6) 



(21) 



where (2n - 1)!! = (2n - l)(2n - 3)(...)(5)(3)(1), r< = 
min(r, n oop ), r> = max(r, r\ oop ), and P™ are the Leg- 
endre polynomials. In practice, we calculate the first 20 
terms of the sum. To prevent the magnetic field from ex- 
tending beyond the surface of the star, A$ is truncated at 
values smaller than 0.5^4^ )ma x, where A0 jlnax is the maxi- 
mum value of A$. Again, the normalization constant Cp 
is chosen to satisfy our choice of /3B,min- 

As a rough guideline to determine whether there is 
sufficient numerical resolution to capture the MRI, we 
compare the grid cell spacing at the position where 
Pb = /3b, mm against the characteristic (minimum un- 
stable) MRI wavelength 



Amri = 



2ttva 

n 



(22) 



where v\ — 2Pb / p is the Alfven velocity of the plasma. 
The ratio of Amri over the grid spacing Aw at £, 1 (Pb = 
/3s,min) is given in Table [l] along with the maximum field 
amplitude inside the star at the initial time. 

2.2.3. Grid, Units, and Parameters 

In all calculations the radial box size is set to three 
times the equatorial radius of the star L ro = 3we- The 
box size along the z-axis is set smaller than the radial di- 
mension L z = \.2we, compensating for the high aspect 
ratio of the stellar profile while also maintaining reason- 
able resolution along the z-axis and accurate multipole 
boundary conditions for the gravity solver. The logarith- 
mic nature of the coordinate system results in a spatial 
resolution of Azn » Q.Q22we and 0.014ra7£ near the ori- 
gin along the radial axis for the 64 3 and 96 3 grids respec- 
tively. This is more than a factor of two improvement in 
resolution over a uniform grid, and increases the radial 
resolution to nearly match the polar (Az w 0.01875ri7E 
and 0.0125-zz7£ for the 64 3 and 96 3 grids). For reference, 
the azimuthal resolution in the equatorial plane of the 
64 3 grid varies from wA(j) w 0.05uje at zcr ~ vje/2 to 
about 10~ 3 vue near the origin. 

Displayed results in all subsequent figures are pre- 
sented in dimensionless code units. We define the equato- 
rial equilibrium surface radius to be the unit coordinate 
dimension (£ = we = 1.7 x 10 8 cm), the dynamical time 
of a spheroid with equatorial surface radius we to be 
the unit interval of time (t = y/w" E /GM w 0.1 s for the 
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r = 5/3 runs), and the maximum central density of the 

star to be the unit density (p — p max ,o = 10 10 g/cm ). 
This in turn implies a scaled magnetic field unit of 
y/GMp maxfi /w E (w 1.7 x 10 14 Gauss for T = 5/3). 

Table [TJ lists all of the numerical simulations we have 
performed in this study, together with their correspond- 
ing parameter sets. The calculations represent a large pa- 
rameter space including the field strength characterized 
by the minimum plasma beta /?B.min, the field orientation 
(toroidal or poloidal) , the equatorial field or current loop 
radius r\ oop , the equation of state polytropic index T, the 
rotational instability parameter 0, the polar to equatorial 
axis ratio A r , and the grid resolution. Also shown in Ta- 
ble [1] is the maximum magnetic field amplitude inside the 
star, and the ratio of the MM wavelength over the local 
radial cell resolution. The prefixes "T" and "P" in our 
naming convention refer to initial toroidal and poloidal 
field orientations, respectively. Runs labeled without ei- 
ther of these prefixes are the baseline unmagnetized cal- 
culations. 



2.3. Gravitational Radiation Diagnostic 

Because the spacetime metric is held fixed in our simu- 
lations, we compute the gravitational radiation produced 
by these systems in the (traceless) quadrupole approxi- 
mation. For an observer located on the symmetry axis, 
the two polarizations of the gravitational wave amplitude 
are written as 



contributions: 



n + = —- 
c 

G2 



4 r V 



c 4 r * xy ' 



(23) 
(24) 



where lij is second time derivative of the reduced 
quadrupole moment of the mass distribution, given in 
flat space Cartesian coordinates x k by 



p(xiXj - -S tJ r 2 ) d 3 x 



(25) 



and r = x 2 + y 2 + z 2 is the spherical distance to the 
source center of mass. The integral is evaluated over the 
entire mass distribution. 

It is well known that the straight-forward procedure 
of computing the components of at each time step, 
and then taking the needed time derivatives of the result 
numerically, produces an unaccept able level of numerica l 
noise in the resulting waveforms ( Finn fc Evans||1990| . 
Instead we differentiate equation (|25| with respect to 
time twice, each time using the evolution equations (|T|) - 
Q to replace the time derivatives that appear in the re- 
sulting integrand. Many of the spatial derivatives which 
are introduced by this procedure can then be eliminated 
using integration by parts. T he resulting expre ssion for 
lij is similar to that given in 
here we include a potential term missing 



(J2000J), except 
trom that ex- 



lij — 



2pViVj - -S ljP V k V k 

2 , 

-pXjdi® - pxidj® + -5 l jpx K d k § 

1 2 

-— (2BiB,j - SijB B k ) - -5ij 



B k B h 
8tt 



+Qij + Qji — T^ijQk 



d 3 3 



(26) 



In the following sections, we present the quantity rh + , to 
scale out the radial dependence, and normalize the wave 
amplitudes by (GM/wec 2 ) 2 so that a direct comparison 
can be made to previous unmagnetized calculations. 

3. RESULTS 

3.1. Toroidal Magnetic Field Configurations 

All calculations with toroidal configurations begin with 
an initial axisymmetric solution in near equilibrium, and 
subsequently evolve in a similar fashion. They display 
the characteristic exponential growth of the even (pri- 
marily m — 2 and 4) non-axisymmetric modes at early 
times, followed by a mode saturation phase during which 
the star develops an elongated bar structure with spi- 
ral arms, then a final bar attenuation phase that redis- 
tributes angular momentum as the star evolves into a 
more axisymmetric configuration again. The early and 
intermediate phases of the bar-mode growth are illus- 
trated in Figure [TJ where we show images of the mass 
density (for the G53Binf case) and magnetic pressure (for 
the TG53B100 case) in the equatorial plane. Both sets 
of images display the same contour levels of the mass 
density (0.5, 0.05, 0.005, 0.0009) normalized to the ini- 
tial maximum mass density p max ,o- The mass densities 
are virtually identical in both magnetized and unmag- 
netized runs. The magnetic pressure is initially confined 
within density contour levels 0.05 and 0.5, a region about 
2cr = 0.25we thick in radius. The magnetic field for this 
configuration is not buoyant but remains mostly confined 
to those level boundaries as the bar mode takes shape. 

Comparing mass density contours from the two simu- 
lations in FigurefT] it is apparent that the introduction of 
magnetic fields (at least for this particular field configu- 
ration) does not affect appreciably the growth of the bar 
mode. This is confirmed in Figure [2] where we compare 
the growth of the amplitude \A m \ for the first few non- 
axisymmetric modes (m = 1, 2, and 4) in the same two 
runs as FigurefT] (G53Binf and TG53B100). For reasons 
of clarity, we do not show the m = 3 curve, but note 
instead that the overall shape and amplitude are very 
similar to the m = 1 mode. The quantity \ A m \ plotted 
in Figure [2] is the same as that introduced in New et al. 
( 2000 1 . It represents the mode amplitude inside a ring 
centered at a fixed cylindrical radius (vd — 0A5tde) in 
the equatorial plane (z = 0) 
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pression, and add artificial viscosity and magnetic field 



where p is the density in the ring, averaged over cells 
along the radial and polar axes as specified by the thick- 
ness of the ring in those directions. 
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TABLE 1 

Simulation Parameters 



S i mil 1 cit ion 


Resolution 


3-a ■ a 
MB, mm 


b 


r (tv) c 


A d 

r 


' loop 




| | max 


G53Binf 


64 3 


OO 


0.295 


0/ o 


ft ^ 


0.208 


— 


— 


- 


G53BinfHR 


96 3 


OO 


0.295 


0/ o 


ft ^ 


0.208 


— 


— 


- 


G20Binf 


64 3 


OO 


0.286 


2 


(1.0) 


0.250 


— 


— 


- 


G30Binf 


64 3 


OO 


0.298 


3 


(0.5) 


0.250 


— 


— 


- 


TG53B1 


64 3 


1 


0.295 


O/o 


^1.0 ) 


0.208 


0.65 


35.0 


5.4 x 10 13 


TG53B10 


64 :i 


10 


0.295 


O/O 


^1.0 ) 


0.208 


0.65 


11.0 


1.7 x 10 13 


TG53B100 


64 3 


100 


0.295 


O/O 


l.O] 


0.208 


0.65 


3.5 


5.4 x 10 12 


TG53B100R3 


64 3 


100 


0.295 


O/O 


l.O] 


0.208 


0.28 


5.6 


2.9 x 10 13 


TG53B100HR 


96 3 


100 


0.292 


O/O 


l.O] 


0.208 


0.65 


5.6 


5.5 x 10 12 


TG53B500 


64 3 


500 


0.295 


5/3 


[1.5) 


0.208 


0.65 


1.6 


2.4 x 10 12 


TG53B500R3 


64 3 


500 


0.295 


5/3 


(1.5) 


0.208 


0.28 


2.7 


1.4 x 10 13 


TG53Ble8 


64 3 


1b8 


0.295 


5/3 


(1.5) 


0.208 


0.65 


0.003 


5.4 x 10 9 


TG20B100 


64 3 


100 


0.286 


2 


(1.0) 


0.250 


0.65 


4.4 


1.8 x 10 13 


TG30B100 


64 3 


100 


0.298 


3 


(0.5) 


0.250 


0.65 


5.2 


4.1 x 10 13 


PG53B10 


64 3 


10 


0.295 


5/3 


(1.5) 


0.208 


0.5 


12.8 


2.3 x 10 13 


PG53B10HR 


96 3 


10 


0.292 


5/3 


(1.5) 


0.208 


0.5 


20.5 


2.7 x 10 13 


PG53B100 


64 3 


100 


0.295 


5/3 


(1.5) 


0.208 


0.5 


4.1 


7.2 x 10 12 


PG53B100HR 


96 3 


100 


0.292 


5/3 


(1.5) 


0.208 


0.5 


6.5 


8.5 x 10 12 


PG53B500 


64 3 


500 


0.295 


5/3 


(1.5) 


0.208 


0.5 


1.8 


3.2 x 10 12 


PG53B500HR 


96 3 


500 


0.292 


5/3 


(1.5) 


0.208 


0.5 


2.9 


3.8 x 10 12 


PG30B100 


64 3 


100 


0.298 


3 


(0.5) 


0.250 


0.65 


6.0 


4.2 x 10 13 



a /3B im in — {P/Pb )mm is the plasma parameter defining the minimum hydro-dynamic to magnetic pressure 
ratio (or maximum magnetic field amplitude). 

(3 — T/\W\ is the rotational instability parameter. 
c r — 1 + 1/N is the adiabatic index for the ideal gas equation of state. 

A r is the polar to equatorial axis ratio of the initial star configuration. 
c r loop is th° radius of the field or current loop in the equatorial plane. 

^ Amri/Act is the ratio of the MRI wavelength to the (cylindrical) radial cell width (Act) at a position 

corresponding to the minimum initial plasma beta (/^B.min) inside the star. 

s I B I max is the maximum magnetic field amplitude inside the star in units of Gauss. 



An obvious result taken from Figure [2] is the strong 
similarity in the early exponential growth and late time 
saturation profiles of both magnetized and unmagnetized 
runs. But there are also several other points of interest in 
Figure [2] For example, the m = 1 mode saturates after 
about five dynamical times to a level that is well below 
those of the to = 2 and 4 modes (about three orders of 
magnitude below m = 2, and two orders below to = 4 at 
their peaks), as is characteristic of the bar mode. Since 
the odd modes tend to encapsulate numerical errors such 
as the center of mass drift and loss of angular momen- 
tum, it is encouraging to see the odd modes saturate in 
time and the characteristic even modes to dominate the 
Fourier spectrum in convincing fashion. We have tracked 
the mass center in our calculations and plot the result 
for a few typical cases in Figure [3] The early drift in the 
center of mass quickly saturates within a few dynamical 
times, and the center of mass thereafter remains mostly 
stationary. The total late-time drift in the position of 
the star mass center is very small, ~ O.Olnj^ for the 64 3 
grids, effectively confined to less than half a cell width in 
the central most highly resolved (smallest zoning) por- 
tion of the grid. The center of mass position is preserved 
even better in the high-resolution 96 3 runs, being con- 
fined to within one quarter of a cell width, ~ 0.0025we- 

Although we do not show the results here, we have 
verified that the modal histories in the magnetized runs 
are very similar to the unmagnetized results also for the 
r = 2 and 3 cases, and both sets look like the T = 5/3 



results in Figure pi The only major difference that we 
have observed ana attributed to the equation of state is 
a systematic shift (or delay) in time required for the even 
bar modes to enter the exponential growth phase. This is 
discussed further in the paragraphs below in the context 
of gravitational wave emissions. 

The ineffectiveness of toroidal magnetic fields is 
demonstrated again by two additional calculations in 
which we varied PB.min to generate different initial field 
amplitudes: a factor of five reduction in run TG53B500, 
and a ten-fold increase in run TG53B10. In both cases, 
toroidal fields do not affect significantly the growth or 
development of the bar mode, in spite of the field am- 
plification observed in Figures [3] and |5] which plot the 
total integrated magnetic energy and the local field am- 
plitude, respectively, as a function of time. The field 
amplitude in Figure [5] is calculated as an azimuthal av- 
erage inside a circular annulus in the equatorial plane 
at radii 0.65we (the initial center of the toroidal field 
loop) and 0.2vue- Note the field amplitude is already 
saturated and does not evolve much at the initial con- 
figuration radius (ri oop = 0.65zue), but grows substan- 
tially at smaller radii. Eventually, the field amplitude 
and energy saturate in time at all radii to roughly the 
same value. This behavior is seen in all of the cases we 
have tried, including the smallest field amplitude case 
(TG53Ble8) in which the field grows from \B\ ~ 3 x 10" 5 
to ~ 0.02 at radius 0.65we by runs end. A confirmation 
of the convergence of these results is provided by the 
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Fig. 1. — Development of the bar mode in the V = 5/3 sim- 
ulations with a magnetic field (TG53B100, right column) and 
without (G53Binf, left column). Rows represent snapshot solu- 
tions at times: t = 12, 14, and 20 in dynamical (code) units 

(t = yJva^,/GM). Contour levels represent mass density in the 

equatorial plane at four values normalized to the initial maximum 
density p max ,o: (0.5, 0.05, 0.005, and 0.0009). Linear color maps 
represent mass density in the case without magnetic fields (left col- 
umn), and magnetic pressure in the case with fields (right). The 
magnetic pressure color scale is fixed for all the images to cover the 
range — * 10 -4 . The mass density color scale is adjusted to the 
maximum density at each time, but does not deviate much from 
— » 1.07 over all displayed sequences. 
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Fig. 2. — Growth of the first few azimuthal Fourier mode am- 
plitudes |j4 m | in the mass density, comparing two T = 5/3 
cases: unmagnetized (G53Binf, dark line types) and magnetized 
(TG53B100, gray lines and symbols). Results are derived by av- 
eraging the density inside a ring centered at w = 0.45toe in the 
equatorial plane, extending one zone deep in both radial and polar 
directions. We do not plot the m = 3 curve here for reasons of 
clarity, but note that it roughly tracks the m = 1 curve in shape 
and amplitude. 
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Fig. 3. — Spherical radial position of the center of mass in a few 
typical unmagnetized runs. The vertical axis is dimensionally nor- 
malized by the initial equatorial surface radius rog. Apart from a 
brief drift at early times, the center of mass movement quickly satu- 
rates and remains confined to a region that is within a half/ quarter 
cell width radius for the 64 3 /96 3 grids. 
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Fig. 4. — Total magnetic energy integrated over the entire grid 
for a few of the T = 5/3 cases. The horizontal axis is plotted in 

dynamical time units ^fu 
sionless code units. 



j 3 e /GM, and the field energy in dimen- 



high- resolution run TG53B100HR which falls between 
the BIO and B500 results, and tracks closely the corre- 
sponding energy and field amplitudes at lower resolution. 
The only difference between the high and lower resolu- 
tion curves is the growth phase in the high resolution 
case is triggered about one dynamical time earlier (ap- 
proximately a 10% shift in time), but saturates at the 
same mean field amplitude. 

The exponential growth of the magnetic energy around 
a time t *~ 11 suggests the onset of a powerful magnetic 
field amplification mechanism. The delayed onset and 
steep exponential growth particularly favor axisymmet- 
ric modes of the MRI. These modes have the shortest 
growth times, but require a poloidal field component to 
act upon. For the initially toroidal configurations, signif- 
icant poloidal fields are not present until the bar mode 
begins redistributing material within the star, explaining 
the delayed onset. Furthermore, it is clear from Figures 
[4] and |5] that the magnetic field saturates to roughly the 
same level at all radii and for all cases, regardless of ini- 
tial amplitudes. The saturation level appears compara- 
ble to the initial field amplitudes in the B500 case: field 
amplitudes in cases initially greater than B500 are dissi- 
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Fig. 5. — Magnetic field amplitude averaged azimuthally inside 
a circular annulus in the equatorial plane between radii 0.65ro£ 
and 0.2-oje- The horizontal axis is plotted in dynamical time units 

A/zOg/GM, and the field amplitude in dimensionless code units of 
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Fig. 6. — Evolution of the mass density weighted average of the 
inverse of the magnetic plasma beta (/3g ) inside the star (satisfied 

by pl Phlox > where p ma x is the maximum density in the 

star at any given time). The high-resolution run TG53B100HR 
(not shown here for reasons of clarity) closely resembles the growth 
curve from the corresponding lower resolution case TG53B100. 

pated through hydrodynamic processes; field amplitudes 
in cases initially less than B500 arc amplified until they 
reach B500 levels before leveling off. Our calculations 
suggest no mechanism exists which can drive field am- 
plitudes above B500 levels, implying that we are safely 
modeling the upper limit of self-generated field strengths. 

Although field amplification does indeed take place in 
all magnetized runs we have performed (medium and 
high resolution), it falls well short of thermal equipar- 
tition so it cannot easily affect the dynamical evolution 
of the star. This is demonstrated in Figure [6] which shows 
the mass density weighted average of the inverse plasma 
beta (1/Pb) inside the star. In all cases the increase of 
magnetic pressure saturates at a level that is significantly 
less than 10% of the thermal pressure averaged across the 
star (with maximum peak values of about 4% for the 64 3 
runs, and 6% for the high resolution 96 3 case). It appears 
that field saturation is determined by universal behavior 
in the partitioning of thermal and magnetic energy, in- 
dependent of initial amplitude. This is true also for the 
r = 2 and 3 cases, both of which result in mean l/fin 
profiles similar to the r — 5/3 results shown in Figure [6j 



Consequently, we do not expect gravitational wave- 
forms to be affected appreciably by toroidal magnetic 
fields, as we demonstrate in Figure [jjfor the r = 5/3 
cases, Figure [8] for r = 2, and FigurelO] for T = 3. Fig- 
ures [7] through [9] plot the quantity rh+ normalized by 
the scale factor {GM/wec 2 ) 2 . Figure m\ also includes 
results from a poloidal initial field configuration (run 
PG30B100) for comparison. Figure [7] (correspondi ng to 
the T = 5/3 cases) closely resembles Figure 9 of New 
et al. ( 2000[ ): we match the amplitude and frequency 
ot oscillations and find that our results are intermediate 
between the "Dl" and "LI" displays in duration and pat- 
tern of the wave signal. It is not until the field strength is 
increased to l/0B,min ~ 1 with locally comparable ther- 
mal and magnetic pressures that we observe amplitude 
deviations of order 30% in Figure [jj 

The first evidence of oscillations in Figure [7] occurs at 
time t ~ 6 which corresponds to the instant when the 
m = 2 mode first begins to dominate the spectral sig- 
nal in Figure [2l The global envelope shape (essentially 
the overall amplitude) of the gravitational wave emission 
tracks nicely the growth and eventual decay of the m = 2 
mode curve in Figure [2] Maximum peaks in both wave 
signals and spectral mode profiles correlate precisely at 
time t ~ 16, and both exhibit comparable rise and decay 
times. Another point of interest in comparing Figures [7] 
-[9] is the apparent trend for the start of the wave signals 
tooe delayed with increasing adiabatic index T (evident 
also in the spectral mode plots). However, we have found 
that the onset of the instability is sensitive to a number of 
numerical factors (e.g., grid resolution, Courant factor), 
and it is difficult to make quantitative conclusions regard- 
ing this effect. For example, the magnetized and unmag- 
netized T = 5/3 high-resolution cases resemble Figure 17] 
but for a slight delay of about 1.5 dynamical times, effec- 
tively a 10% temporal shift in the waveform. However, 
other aspects of the waveforms are similar between the 
higher and lower resolution cases: the magnetized results 
are essentially identical to the unmagnetized waveforms, 
and the wave amplitudes agree nicely. 

Even though the gravitational wave amplitudes are 
fairly consistent with no obvious correlation with T (ap- 
proximately 0.45, 0.35, and 0.44 for T = 5/3, 2 and 3, 
respectively), the wavelength of perturbations between 
the two biggest wave crests and the burst duration (be- 
tween leading and trailing wave crests with amplitude 
larger than 0.1) do appear to increase monotonically with 
r. In particular, we find wavelengths of approximately 
3.lt, 4.0t, and 4.8? for T = 5/3, 2, and 3, a fractional in- 
crease (in wavelength) of about 20% between T = 2 and 
r = 5/3, and about 55% between T = 3 and T = 5/3. 
Trends in both a mplitudes and wavelengths are consis- 
tent with those of |Houser fc Centrella (1996), who find 
the amplitude is independent ot the poiytropic index, but 
r = 2 (r = 3) fluids generate gravitational waves with 
20% (58%) longer wavelengths than T = 5/3. 

We cannot quantify the burst duration as easily since 
the late time solutions can be inaccurate due to buildup 
of numerical errors, and sensitivity to computational pa- 
rameters (e.g., Courant factor, artificial viscosity con- 
stants). However, a comparison of Figures [jj - [9] clearly 
shows a lengthening of the pulse duration by more than 
50% as the m odels grow stiffer . Such behavior was ob- 
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Fig. 7. — Gravitational waveforms (rh+) extracted from several 
magnetized and unmagnetized T = 5/3 cases. The horizontal time 

axis is displayed in dynamical code units <J vj^/GM, and the wave 

amplitude is normalized by (GM/s7gc 2 ) 2 . 
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Fig. 8. — Gravitational waveforms with magnetic fields (dashed 
line) and without (solid line) for T = 2. 



& Tohline ( 1988 1 in their unmagnetized studies: stiffcr 
polytropes produce more elongated bars, rotate more 
slowly, and undergo more periods of spiral arm ejection 
and core recontraction, resulting in longer bursts of grav- 
itational wave signals. The addition of toroidal magnetic 
fields to the stellar profile does not appear to affect these 
behaviors appreciably as the burst duration is generally 
shorter than the timescale ts for magnetic braking to 
take effect: ts ~ zue/va > 25-50 dynamical (code) units 
at the field saturation time for all of the cases we have 
considered. 

3.2. Poloidal Magnetic Field Configurations 

Unlike the cases which start from initially toroidal 
magnetic field configurations, simulations performed of 
models which begin with poloidal fields evolve differently 
than unmagnetized cases if the initial field amplitude is 
large. In particular, we found little qualitative difference 
between all the toroidal calculations at both low (64 3 ) 
and high (96 3 ) grid resolutions, regardless of initial field 
amplitude . Fo r this reason, all of the results presented 
in Section [3 . 1 1 corresponded to low-resolution cases, with 
the added advantage of allowing for many different pa- 
rameter combinations to be investigated. However, for 
the poloidal cases, we observed an increased sensitivity to 
spatial resolution as well as a growing impact on bar for- 
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Fig. 9. — Gravitational waveforms with magnetic fields (dashed 
lines) and without (solid line) for T = 3. 

mation with increasing field strength, especially for the 
two largest amplitude cases (PG53B100 and PG53B10). 
All poloidal results presented in this section are therefore 
shown at high (96 3 ) grid resolution. 

To illustrate the effect of the magnetic field on bar 
formation, Figure [10] shows images of the mass den- 
sity (for the G53BRifHR case) and magnetic pressure 
(for the cases PG53B100HR with /3 B ,min = 100, and 
PG53B10HR with /3b, mm = 10), all in the equatorial 
plane. All three sets of images display the same con- 
tour levels of the mass density (0.5, 0.05, 0.005, and 
0.0009), normalized to the initial maximum mass den- 
sity pmax,o- Images making up the left column for the 
case without a magnetic field are shown at times t =15, 
17, and 21; images in the center column are shown for 
case PG53B100HR at times t =17, 19, and 23; and im- 
ages in the right column are shown for case PG53B10HR 
at times i=15, 17, and 19, all in dynamical (code) units. 
The PG53B100HR images are shown at later times than 
those for the unmagnetized case to illustrate that a bar 
forms with a shape similar to the no-field case, but de- 
layed approximately two dynamical times. For the higher 
magnetic field strength run, PG53B10HR, there is no in- 
dication that a structure resembling a bar is going to 
form at this resolution, a point we return to below. 

The effect of an initially poloidal magnetic field on the 
growth of the bar mode is quantified in the nonaxisym- 
metric am plitu des \A m \, c omp ute d as described above 
in Section |3.1[ In Figures [TT] and [l2j we plot the evo- 
lution of the amplitudes \A m \ of the m = 2 and m = 4 
modes, respectively, for runs G53BInfHR, PG53B500HR, 
PG53B100HR, and PG53B10HR. It is apparent that for 
the lowest field case (PG53B500HR), the magnetic field 
has little effect on the growth of either mode, but increas- 
ing the field amplitude systematically suppresses the dif- 
ferent modes. For example, the intermediate field case 
(PG53B100HR) exhibits growth similar to the no-field 
case but with a temporal delay and slightly smaller mode 
amplitudes, while modes in the highest field amplitude 
case (PG53B10HR) are suppressed by approximately two 
orders of magnitude relative to the no-field case. Al- 
though not shown, we note that the m = 1 and m = 3 
modes are similar to those observed in the toroidal cases, 
in that they saturate at levels which are well below the 
characteristic m = 2 and 4 profiles of the bar mode. 

We expect the gravitational waveforms produced in 
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Fig. 10. — Development of the bar-mode in T = 5/3 simulations 
with initially poloidal magnetic fields (PG53B100HR, center col- 
umn, and PG53B10HR, right column) and without (G53BInfHR, 
left column). Rows represent snapshot solutions at different times: 
t = 15, 17, and 21 for the case without magnetic fields, t = 17, 19, 
and 23 for the PG53B100HR case, and t = 15, 17, and 19 for the 

PG53B10HR case, all in dynamical (code) units (i = ^Jm 3 E /GM). 
Contour levels represent mass density in the equatorial plane at 
four values normalized to the initial maximum density p m ax,o ; (0-5, 
0.05, 0.005, and 0.0009). Linear color maps represent mass density 
in the case without magnetic fields (left column), and magnetic 
pressure in the cases with fields (center and right columns). The 
magnetic pressure color scale is fixed for all the images to cover the 
range — » 10 . 
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Fig. 11. — Growth of the m = 2 azimuthal Fourier mode am- 
plitudes \A m \ in the mass density, comparing four T = 5/3 cases: 
unmagnetized (G53BInfHR, dark line) and magnetized, with ini- 
tially poloidal magnetic fields (PG53B500HR, long dashed line; 
PG53B100HR, short dashed line; and PG53B10HR, dotted line). 
Results are derived by averaging the density inside a ring centered 
at zu = OAbvjE in the equatorial plane and extending one zone 
deep in both radial and polar directions. 

each of these runs to reflect what is seen in the den- 
sity profiles and the mode amplitude plots. In Figure 13 
we show the quantity rh + normalized by the scale fac- 
tor (GM/w E c 2 ) 2 , for runs G53BInfHR, PG53B500HR, 
PG53B100HR, and PG53B10HR. For the plot shown, the 
low-amplitude field case PG53B500HR has been shifted 
back 0.2 dynamical times, and the moderate-amplitude 
field case PG53B100HR has been shifted 1.4 dynami- 
cal times so that the maximum of the first large peaks 
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Fig. 12. — Growth of the m = 4 azimuthal Fourier mode am- 
plitudes \A m \ in the mass density, comparing four T = 5/3 cases: 
unmagnetized (G53BInfHR, dark line) and magnetized, with ini- 
tially poloidal magnetic fields (PG53B500HR, long dashed line; 
PG53B100HR, short dashed line; and PG53B10HR, dotted line). 
Results are derived by averaging the density inside a ring centered 
at zu — 0.45rog in the equatorial plane and extending one zone 
deep in both radial and polar directions. 
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Fig. 13. — Gravitational waveforms (r/i+) extracted from four 
r = 5/3 cases: unmagnetized (G53BInfHR, dark line) and magne- 
tized, with initially poloidal magnetic fields (PG53B500HR, long 
dashed line; PG53B100HR, short dashed line; and PG53B10HR, 
dotted line). In this plot, the waveforms produced in the 
PG53B500HR and PG53B100HR simulations have been shifted in 
time so that the peaks in the waveforms line up with those of the 
unmagnetized case. 

line up. As expected, waveforms exhibit systematically 
smaller amplitudes and greater temporal delays with in- 
creasing initial field amplitude. 

In Figure |14| we plot the inverse plasma beta 
(1//3b), azimuthally averaged over the ring at w = 
0.65h7e in the equatorial plane, for runs PG53B500HR, 
PG53B100HR, and PG53B10HR. For the lowest field 
case (PG53B500HR), in which the dynamics are affected 
only slightly, the magnetic pressure is never more than 
about 4% of the fluid pressure in this ring. However, for 
the highest field case (PG53B10HR), in which the bar 
mode is completely suppressed, the magnetic pressure is 
already just below 4% of the fluid pressure in the initial 
configuration, and peaks at around 80% of the fluid pres- 
sure within five dynamical times before leveling off again 
at ~4%. Under these conditions, we would expect the 
magnetic field to have a significant effect on the overall 
evolution of the star. 

We can understand the immediate onset of the field 
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Fig. 14. — Inverse of the magnetic plasma beta (1//3b) az- 
imuthally averaged over a ring in the equatorial plane at the radius 
0.65h7e. 

amplification observed in Figure [14] by considering the 
evolution of i ndividual components of the magnetic field. 
In Figure [15] we plot 8tt times the magnetic pressure, as 
well as the contributions to this quantity from each of the 
three cylindrical components of the magnetic field, all az- 
imuthally averaged over the ring at vo = 0.65we in the 
equatorial plane, for runs PG53B500HR, PG53B100HR, 
and PG53B10HR. For all three cases, we see that the in- 
crease in the overall magnetic pressure tracks the increase 
in the azimuthal contribution to that pressure. This be- 
havior, in addition to the immediate field amplification, 
is characteristic of the f2-dynamo, which one would ex- 
pect to be active with a poloidal magnetic field present. 
This effect can be seen qualitatively in Figure 



16 



which we show magnetic field lines at early times for the 
PG53B100HR case. The magnetic field lines, as well as 
mass density contours, are shown at times t — 0, 1, 2, and 
3t. The field lines show immediate and dramatic stretch- 
ing in the azimuthal direction, c ontr ibut ing to the field 



amplification observed in Figures [14| and 15 

Another important point from Figure 14 is that the 
late time saturation level of 1/Pb is approximately the 
same in all three simulations. A similar behavior and 
level of saturation are evident in Figure [6] for the toroidal 
cases. This remarkable independence from the initial 
field strength and orientation is important, as it indi- 
cates that a common and robust physical mechanism is 
at play setting the final saturation level. Also, we see 
from Figure [15] that the final magnitude of the azimuthal 
component (g^B^B^) is approximately the same in all 
three simulations, and the second most significant com- 
ponent, the radial one (g nn B v B v ), always finishes within 
an order of magnitude of the azimuthal component. This 
suggests that these two components have either achieved 
a static final configuration common to all the simulations, 
or they arc feeding off of one another and exchanging 
energy with the gas through some equilibrium process. 
The latter is consistent with expectations of the MRI, 
which gradually takes over as the fi-dynamo saturates, 
and which is expected to be active in all of our stellar 
models according to criterion (Equation (fl8|)). 

We end this section with a few words of caution regard- 
ing the convergence of results for those cases run with 
poloidal initial fields. The smallest B500 case is reason- 
ably converged between low (64 3 ) and high (96 3 ) grid res- 
olutions, resulting in nearly identical gravitational wave 
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Fig. 15. — We show Sit times the magnetic pressure, as well as 
the contributions to this quantity from each of the three cylindrical 
components of the magnetic field, all azimuthally averaged over the 
ring at = 0.65ote in the equatorial plane, for runs PG53B500HR 
(top), PG53B100HR (center), and PG53B10HR (bottom). 

frequencies and maximum peak differences of about 25% 
in amplitude. However, the intermediate B100 case ex- 
hibited vastly different solutions at low and high reso- 
lutions. The bar mode was essentially suppressed com- 
pletely at low resolution, but managed to form nicely 
at 96 3 zones, demonstrating the sensitive and demand- 
ing nature of these calculations. We cannot therefore be 
certain that the largest amplitude (BIO) case is fully con- 
verged and does not form a bar. It would require even 
larger computational grids to confirm this result, some- 
thing that is beyond our current allocation resources. At 
issue is the resolution of the winding of the field lines 
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Fig. 16.— Magnetic field lines for the PG53B100HR simulation, 
at times and 1 t (top row), and 2 and 3 t (bottom row). Also 
shown are mass density contours at three values, normalized to the 
initial maximum density p max ,o: (0.5, 0.05, 0.005). 

as the star rotates, the accurate capturing of reconnec- 
tion events, and mostly the severe restrictions placed 
on the stability time step from magnetically dominated 
MHD and fast Alfven velocities as the field works its way 
through the stellar atmosphere. 

4. CONCLUSIONS 

We have studied the growth of the dynamical bar-mode 
instability in differentially rotating magnetized neutron 
stars through a set of numerical Newtonian MHD cal- 
culations. The calculations explored both toroidal and 
poloidal initial field distributions of differing strengths, 
as well as the role of the equation of state. 

For our initially toroidal field configurations, field am- 
plification always saturated at a level insufficient to 
strongly affect the dynamics of the bar mode. Even the 
most extreme case (TG53B1, /3B,min = 1), with equal 
initial magnetic and hydrodynamic pressures within the 
toroidal loop, gave only a ~ 30% change in gravitational 
wave properties. Otherwise evolution proceeded quite 
similarly to our unmagnetized reference cases. 

The effects were larger for the poloidal configurations. 
In the most extreme case considered (model PG53B10HR 
with /3b. mm = 10), the magnetic field was sufficient to 
completely suppress the formation of a bar. However, 
in that case, we started with the magnetic field already 
within a factor of ~30 of being dynamically more im- 
portant than thermal pressure in some portions of the 



star. It then only took a few dynamical times of field 
growth for /3b to approach unity, a timescale that is much 
shorter than the bar deformation time for azimuthal 
Fourier modes to reach appreciable amplitudes. Also, 
due to the computationally demanding nature of evolv- 
ing strong poloidal fields, we cannot verify that this re- 
sult will not change with increased grid resolution. For 
less extreme initial conditions, the effects of including 
poloidal field components were consistently more mod- 
est with decreasing initial field amplitude. These lower 
amplitude calculations are also less demanding compu- 
tationally and can, like the toroidal cases, be more easily 
checked for convergence. The principle effects of intro- 
ducing poloidal fields are in systematically delaying the 
onset of the bar mode and in suppressing both Fourier 
mode and gravitational wave amplitudes, all of which 
become negligible with initial field amplitudes below the 
threshold of l//3 B ,mm < 10~ 2 . 

Overall, our results suggest that the effect of magnetic 
fields on the emergence of the bar-mode instability in 
neutron stars is not likely to be very significant. Particu- 
larly considering that collapse progenitor models predict 
realistic field configurations that are dominantly toroidal 
in nature with toroidal and poloi dal compone nts of order 
10 10 G and 10 6 G, respectively ( |Heger et al.]|2005| , sub- 
stantially below most of the field configurations we have 
considered. Thus, except in special cases where neutron 
stars are born very highly magnetized, we might still ex- 
pect them to be good gravitational wave sources if their 
rotational kinetic energies exceed the critical bar-mode 
instability parameter. 



Computations were performed at the Barcelona Su- 
percomputing Center (BSC) under activity AECT-2007- 
3-0002, the Lawrence Livermore National Laboratory 
(LLNL), the College of Charleston (CoC), and at the 
High Performance Academic Computing Environment 
at Washburn University (HiPACE). This work was per- 
formed in part under the auspices of the U.S. Depart- 
ment of Energy by Lawrence Livermore National Labo- 
ratory under contract no. DE-AC52-07NA27344. P.C.F. 
gratefully acknowledges the support of the College of 
Charleston 4th Century Initiative and the South Car- 
olina Space Grant Consortium. J.A.F. acknowledges fi- 
nancial support from the Spanish Ministry of Education 
and Science (AYA 2007-67626-C03-01). 



REFERENCES 



Anninos, P., & Fragile, P. C. 2003, ApJS, 144, 243 
Anninos, P., Fragile, P. O, & Murray, S. D. 2003, ApJS, 147, 177 
Anninos, P., Fragile, P. O, & Salmonson, J. D. 2005, ApJ, 635, 
723 

Baiotti, L., Giacomazzo, B., & Rezzolla, L. 2008, Phys. Rev. D, 
78, 084033 

Baiotti, L., Pietri, R. D., Manca, G. M., & Rezzolla, L. 2007, 

Phys. Rev. D, 75, 044023 
Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 
Burrows, A., Dessart, L., Livne, E., Ott, C. D., & Murphy, J. 

2007, ApJ, 664, 416 
Cerda-Duran, P., Font, J. A., Anton, L., &; Miiller, E. 2008, A&A, 

492, 937 

Chandrasekhar, S. 1969, Ellipsoidal figures of equilibrium (The 
Silliman Foundation Lectures; New Haven, CT: Yale Univ. 
Press) 



Dedner, A., Kemm, F., Kroner, D., Munz, T., C. D. Schnitzer, & 
Wesenberg, M. 2002, J. Comput. Phys., 175, 645 

Dimmelmeier, H., Ott, C. D., Marek, A., & Janka, H.-T. 2008, 
Phys. Rev. D, 78, 064056 

Finn, L. S., & Evans, C. 1990, ApJ, 351, 588 

Fragile, P. C, Anninos, P., Gustafson, K., & Murray, S. D. 2005, 

ApJ, 619, 327 
Hachisu, I. 1986, ApJS, 61, 479 

Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 
Houser, J. L., & Centrella, J. M. 1996, Phys. Rev. D, 54, 7278 
Houser, J. L., Centrella, J. M., & Smith, S. C. 1994, Phys. Rev. 
Lett., 72, 1314 

Jackson, J. D. 1975, in Classical Electrodynamics (2nd ed.; New 

York: Wiley) 
Liu, Y. T. 2002, Phys. Rev. D, 65, 124003 



Bar-mode instability in magnetized neutron stars 



13 



New, K. C. B., Centrella, J. M., & Tohline, J. E. 2000, 

Phys. Rev. D, 62, 064019 
Ott, C. D., Burrows, A., Thompson, T. A., Livne, E., & Walder, 

R. 2006, ApJS, 164, 130 
Ott, C. D., Dimmelmeier, H., Marck, A., Janka, H.-T., Hawke, I., 

Zink, B., & Schnetter, E. 2007, Phys. Rev. Lett., 98, 261101 
Shibata, M., Baumgarte, T. W., & Shapiro, S. L. 2000, ApJ, 542, 

453 



Shibata, M., Karino, S., & Eriguchi, Y. 2002, MNRAS, 334, L27 
Shibata, M., Karino, S., & Eriguchi, Y. 2003, MNRAS, 343, 619 
Spruit, H. C, & Phinney, E. S. 1998, Nature, 393, 139 
Williams, H., & Tohline, J. 1988, ApJ, 334, 449 
Woosley, S. E., & Heger, A. 2006, ApJ, 637, 914 



